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Abstract: We compute the relative-order-v 4 contribution to gluon fragmentation into 
quarkonium in the 3 5i color-singlet channel, using the nonrelativistic QCD (NRQCD) 
factorization approach. The QCD fragmentation process contains infrared divergences 
that produce single and double poles in e in 4 — 2e dimensions. We devise subtractions that 
isolate the pole contributions, which ultimately are absorbed into long-distance NRQCD 
matrix elements in the NRQCD matching procedure. The matching procedure involves 
two-loop renormalizations of the NRQCD operators. The subtractions are integrated over 
the phase space analytically in 4 — 2e dimensions, and the remainder is integrated over the 
phase-space numerically. We find that the order-v 4 contribution is enhanced relative to 
the order-u° contribution. However, the order-f 4 contribution is not important numerically 
at the current level of precision of quarkonium-hadroproduction phenomenology. We also 
estimate the contribution to hadroproduction from gluon fragmentation into quarkonium 
in the 3 Pj color-octet channel and find that it is significant in comparison to the complete 
next-to-leading-order-in-a s contribution in that channel. 
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1 Introduction 

Heavy- quarkonium production in hard-scattering collisions has a long and rich history of 
experimental measurements and theoretical calculations [1]. Intense efforts in this area 
are expected to continue as the Large Hadron Collider (LHC) makes available data with 
unprecedented momentum transfers and statistics. 

In recent years, a great deal of theoretical effort has been focused on the nonrelativistic 
QCD (NRQCD) factorization approach [2] to calculations of quarkonium production rates. 
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In this approach, it is conjectured that the inclusive quarkonium production cross section 
at large transverse momentum (pj>) can be written as a sum of products of short-distance 
coefficients and long-distance matrix elements (LDMEs): 

a(IT) = 2>n(MA)<O|0j?(/iA)|O>. (l.l) 

n 

Here, fi\ is the factorization scale, which is the cutoff of the effective field theory NRQCD. 
A short-distance coefficient F n (fj,\) is, in essence, the partonic cross section to produce 
a heavy-quark-antiquark (QQ) pair with certain quantum numbers, convolved with par- 
ton distribution functions. The short-distance coefficients can be calculated as pertur- 
bation series in the strong-coupling constant a s . A production LDME (O|0 n (/^y)|0) is 
the probability for a QQ pair with certain quantum numbers to evolve into a particular 
heavy-quarkonium state. It is expressed as the vacuum expectation value of a four-fermion 
operator 

0% (Ma) = (OlxW^CP^^xlO), (1-2) 

where ift and x are two-component (Pauli) fields that create a heavy quark and a heavy 
antiquark, respectively, and K n and n' n are combinations of Pauli and color matrices. 1 

V H{P) = £ \H{P) + X)(H(P) + X\ (1.3) 
x 

is a projection onto a state consisting of a quarkonium H, with four-momentum P, plus 
anything. VniP) contains a sum over any quarkonium polarization quantum numbers that 
are not specified explicitly. The NRQCD LDMEs are evaluated in the rest frame of the 
quarkonium, in which P = (M, 0), where M is the quarkonium mass. In the remainder of 
this paper, we suppress the momentum argument of Vjj^ in NRQCD LDMEs. 

The production LDMEs for the evolution of color-singlet QQ pairs into a quarkonium 
state are related to the color-singlet quarkonium decay LDMEs. These color-singlet pro- 
duction LDMEs can be determined from comparison of theory with quarkonium production 
or decay data or from lattice QCD calculations. However, the production LDMEs for the 
evolution of color-octet QQ pairs into a quarkonium state can be determined, at least 
at present, only through comparison of theory with experimental quarkonium- production 
data. 

Complete calculations of short-distance coefficients in the NRQCD factorization ap- 
proach now exist through next-to- leading order (NLO) in a s for production of the J/tp 
and the ip(2S) in e + e~ collisions, in ep collisions, and in pp and pp collisions [1, 5-15]. 
These calculations include the contributions from all of the color-octet channels through 
relative order t> 4 , as well as the contribution of the color-singlet channel at leading order 
(LO) in v. Specifically, the calculations include the contributions of the S\, 1 Sq and 3 Pj 
color-octet channels and the 3 Si color-singlet channel at the leading nontrivial order in v 
in each channel. Here, v is half the relative velocity of the heavy quark and the heavy 

1 It was pointed out by Nayak, Qiu, and Sterman that gauge invariance requires that the definitions of 
the NRQCD LDMEs include Wilson lines that run from the quark and antiquark fields to infinity [3, 4]. 
For simplicity, we have omitted these Wilson lines here. 
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antiquark in the quarkonium rest frame, v 2 ~ 0.22 for the J/ip, and v 2 ~ 0.1 for the 
T. These theoretical results are generally compatible with experimental measurements of 
quarkonium production cross sections. However, significant discrepancies remain between 
theoretical predictions for quarkonium polarization and experimental measurements [13- 
15]. These discrepancies might point to as-yet-uncalculated theoretical contributions, to 
experimental difficulties, to a failure of convergence of the NRQCD series in a s or v, or to 
a failure of the NRQCD factorization conjecture itself. 

The calculations at NLO in a s have revealed very large corrections in that order to 
the 3 -Pj color-octet channel and the 3 Si color-singlet channel. Large corrections have also 
been found in a calculation of the real-emission contributions to T hadroproduction at 
next-to- next-to-leading order in a s [16]. The large corrections are the result of kinematic 
enhancements of the higher-order cross sections at high p T relative to the LO cross sections. 
The size of these corrections has cast some doubt on the convergence of the perturbation 
series. 

It has been suggested recently that the large higher-order corrections to quarkonium 
production can be brought under control by re-organizing the perturbation series according 
to the Pt behavior of the various contributions [17]. In this approach, the cross section can 
be shown to factorize into convolutions of hard-scattering cross sections with fragmentation 
functions. The factorization holds up to corrections of relative order rri^/p 4 -,, where m c is 
the charm-quark mass. The factorized cross section consists of a leading contribution, 
which arises from single-particle fragmentation into a quarkonium and falls as l/p^ in the 
partonic cross section and a first subleading contribution, which arises from two-particle 
fragmentation into a quarkonium and falls as ~Vjp\ in the partonic cross section. If NRQCD 
factorization holds, then the various fragmentation functions can be expressed in terms of a 
sum of products of short-distance coefficients and NRQCD LDMEs. This picture has been 
shown to account for the large corrections at NLO in a s in the 3 Si color-singlet channel 
[18]. 

In this paper, we compute the NRQCD short-distance coefficient for gluon fragmenta- 
tion to a 3 Si color-singlet QQ pair in relative order v 4 . The short-distance coefficients for 
gluon fragmentation in this channel have been computed in relative order v° [19, 20] and 
relative order v 2 [21]. In both cases, the contributions are not important phenomenologi- 
cally. Nevertheless, it is worthwhile to consider the order-v 4 contribution for two reasons. 
First, this contribution is interesting theoretically because it is in order v A that the 3 Si 
color-singlet fragmentation channel first develops soft divergences in full QCD. As we shall 
see, these soft divergences in full QCD correspond to soft divergences in the LDMEs for the 
3 Si and 3 Pj color-octet channels and cancel in the short-distance coefficients, as is required 
by NRQCD factorization. A second motivation for examining the order-u 4 contribution is 
that it is potentially large. Contributions from gluon fragmentation into the 3 Si and 3 Pj 
color-octet channels are known to be important phenomenologically. The 3 Si color-singlet 
channel mixes with these channels in order v 4 , and the partitioning of the various contri- 
butions is controlled by single and double logarithms of the factorization scale. Therefore, 
it is plausible that the order-t; 4 contributions to the color-singlet channel could be large. 

Our method of calculation is based on the Collins-Soper definition [22] of the frag- 
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mentation function for a gluon fragmenting into a quarkonium. We assume that NRQCD 
factorization holds, that is, that the fragmentation function can be decomposed into a 
sum of products of short-distance coefficients and NRQCD LDMEs. We then compute 
the full-QCD fragmentation functions for a gluon fragmenting into free QQ states with 
various quantum numbers. The ultimate aim is to match these full-QCD fragmentation 
functions to the corresponding NRQCD fragmentation functions in order to determine the 
NRQCD short-distance coefficients. Some additional details of this approach can be found 
in Ref. [21]. 

This paper is organized as follows. We give the Collins-Soper definition of the frag- 
mentation function in Sec. 2. Section 3 contains the NRQCD factorization formula for 
the fragmentation function and also contains a discussion of the NRQCD LDMEs and 
short-distance coefficients that are relevant through relative order v 4 . The kinematics and 
variables that we use in our calculation are described in Sec. 4. In Sec. 5, we discuss the 
calculation of the fragmentation processes in full QCD. As we have mentioned, an impor- 
tant feature of the present calculation is that soft divergences arise in the 3 S± color-singlet 
channel in both full QCD and NRQCD. These divergences ultimately cancel in the short- 
distance coefficients when we carry out the matching between full QCD and NRQCD. In 
both full QCD and NRQCD, we regulate the divergences dimensionally. In the case of 
full QCD, we devise subtractions that remove the divergent terms from the integrand, and 
we compute the subtraction contributions analytically. This computation is described in 
Sec. 5.3. After we remove the subtraction terms, we can calculate the remainder of the 
full-QCD contribution in four dimensions, carrying out the integration numerically. We 
compute the relevant NRQCD LDMEs for free QQ states analytically in dimensional reg- 
ularization. These calculations are described in Sec. 6. We also determine the evolution 
equations for the LDMEs and find a discrepancy with a result in Ref. [23]. In Sec. 7, 
we match the NRQCD and full-QCD fragmentation functions to obtain the short-distance 
coefficients, and we present numerical results for them in Sec. 8. Finally, in Sec. 9, we 
summarize our results. 

2 Collins-Soper definition of the fragmentation function 

Here, and throughout this paper, we use the following light-cone coordinates for a four- 
vector V: 



V 



(v + ,v-,v ± ) = (v+,v-,v\v 2 ) 



(2.1a) 
(2.1b) 
(2.1c) 




The scalar product of two four- vectors V and W is then 



v -w = v + w~ + v~w 



■+ 



v± ■ w±. 



(2.2) 
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The Collins-Soper definition for the fragmentation function for a gluon fragmenting into a 
hadron (quarkonium) H [22] is 

„<2-3 /"+oo 

D[5 - = 2^(4" m- 2) L ^~ e ~" *~ 

x (0)^+^(0)^(0, 0, ± ) cb V H{P) £(0, x-,0 ± ) ba G+ u (0, x~,0 ± )|0) . 

(2.3) 

Here, 2 is the fraction of the gluon's + component of momentum that is carried by the 
hadron, G^ u is the gluon field-strength operator, k is the momentum of the field-strength 
operator, /ua is the factorization scale, and d = 4 — 2e is the number of space-time dimen- 
sions. There is an implicit average over the color and polarization states of the initial gluon. 
The projection Vh{p) is given in Eq. (1.3). The fragmentation function is evaluated in the 
frame in which the hadron has zero transverse momentum: P = [zk + , M 2 / (2zk + ), 0j_]. 
The operator £(0,x _ ,0j_) is a path-ordered exponential of the gluon field: 



£(0,x ,0±)ba = Pexp 



/•oo 

+ig s I dz-A + (0,z-,0 ± ) 

J x~ 



(2.4) 



ba 



where g s = y/4iTa s is the QCD coupling constant and A fl (x) is the gluon field. Both and 
G^ u are SU(3) matrices in the adjoint representation. The expression (2.3) is manifestly 
gauge invariant. We use the Feynman gauge in our calculation. 

The Feynman rules for the perturbative expansion of Eq. (2.3) are given in Ref. [22]. 
The quantity £(0, x~,0±) appears in the Feynman rules as an eikonal line. Owing to the 
charge-conjugation properties of the QQ states that we consider and the Landau- Yang 
theorem [24, 25], gluon attachments to the eikonal lines from £ (0, x~ , 0j_) do not appear in 
our calculation. Hence, we need only the standard QCD Feynman rules, an overall factor 

z d— 3^+ 

Qrag = 27r(JV|-l)(d-2) (2 " 5) 

from Eq. (2.3), and the special Feynman rule for the vertex that creates a gluon and an 
eikonal line. That vertex is shown in Fig. 1. Its Feynman rule, in momentum space, is a 
factor 

9 ua ~ ^f) ^b, (2.6) 

where k is the sum of the momenta of the gluon and the eikonal line, Q is the momentum 
of the gluon, a is the polarization index of the gluon, and a and b are the color indices, 
respectively, of the gluon and the eikonal line. In the absence of interactions with the 
eikonal lines, k = Q. n is a light-like vector whose components are given by n = (0, 1, 0j_). 

The final-state phase space that is implied by Eq. (2.3) is 
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a, a 




Figure 1. Feynman diagram for the vertex that creates a gluon and an eikonal line. The circle 
represents the operator G^ v , which creates a gluon with momentum Q and polarization and color 
indices a and a, respectively, b is the color index for the eikonal line. The operator momentum k 
is the sum of Q and the momentum of the eikonal line. 

where S is the statistical factor for identical particles in the final state, fcj is the momentum 
of the ith final-state particle, and the product is over all of the final-state particles except 
H. We use nonrelativistic normalization for the state H, and so a factor 2M appears in 
the phase space in order to cancel the relativistic normalization of H in the definition (2.3). 
We use relativistic normalization for all particles other than H. 

3 NRQCD factorization 

We assume that the fragmentation function for a gluon fragmenting into a quarkonium H 
satisfies NRQCD factorization. Then, in analogy with Eq. (1.1), we have 

D[g^H](z) = J2dn(z)(0\O^\0), (3.1) 

n 

where the (0|O^|0) are NRQCD LDMEs and the d n(z) are the fragmentation short-distance 
coefficients. We have suppressed the dependences of d n (z) and (O|0^|O) on the factorization 
scale fiA- In discussing specific cases, we use the notation (O|0^ n /( 2s+1 ^ c ')|0) for the 
LDMEs and the notation d nv ^ n i( 2s+1 l^){z) for the short-distance coefficients, where 2s+1 lj 
is the standard spectroscopic notation for the angular-momentum quantum numbers of the 
corresponding NRQCD operator, c is the color quantum number of the NRQCD operator 
(1 or 8). n v is the order in v, relative to the leading order, of the field operators and 
derivatives in n ,, excluding factors of v from the projection onto the final state H. n' 
is an integer that is used to distinguish operators that have the same quantum numbers 
and order in v. We denote the contribution of order v k to D[g — >• H] by Dk[g — > H]. 

We can write the fragmentation functions for gluon fragmentation into free QQ states 

as 

D[g -> QQ\{z) = ^d n (z)(Q\O%Q\0). (3.2) 

n 

Since the short-distance coefficients d n (z) are independent of the specifics of the hadronic 
states, the d n (z) in Eq. (3.2) are identical to the d n (z) in Eq. (3.1). We determine the d n (z) 
by computing the left side of Eq. (3.2) in full QCD and comparing it with the right side, in 
which the free QQ LDMEs are computed in NRQCD. Since we choose a factorization scale 
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//A of order the heavy-quark mass m, we can carry out this computation in perturbation 
theory. We will denote the contribution of order v k to D[g — > QQ] by Dk[g — > QQ]- 

If H is a 3 Si quarkonium state, such as the J/tp, then, in LO in v, we must consider 
the LDME 

(o\o^es?)\o) = (0|xVV v u </>V X |0>. (3.3) 
In relative order t> 2 , we must consider the LDME 

(0|Of ( 3 ^ 1] )|0) = 1<0| X V(_!£) 2 ^ ^V x + H . c .|0>. (3.4) 

In relative order u 3 , we must consider the LDME 

(oi^C^f ] )|0) = (o|x f r> ^r a x |o). (3.5) 

In relative order u 4 , we must consider the LDMEs 

(0\O^esP)\0) = (0| X V(_|£)^ p H v ,t^ ( _^ ) 2 x | 0)) (3 . 6a) 

(0|O|f 2 ( 3 5l 1] )|0) = I(0| x V(-iS)V ^ V>V X + H. c.|0), (3.6b) 

(o|o|f3( 3 5[ 1] )|o) = ^o| x V^^vV(S-^ + ^-S)x 



-x 



^^{D-gE + gE-b^TH ^V X |0), (3.6c) 

(O|0,f ( 3 sf ] )|0> = (0| X Vr> V H ^V^xIO), (3.6d) 

(O|0 O *( 3 P O [1] )|O> = JL-(0| x t(_«3 p H ^t(_|3 . CT)x | ), (3.6e) 

<0|O H ( 3 pj 1] )|0> = (0| x t(_|B[Vi])^ Ph ^(-0^ i] )x\O), (3.6f) 

<O|0 o H ( 3 pj 1] )|O> = (0| x t(_|B(V'))^ Vu ^(-!2M)x|0>. (3.6g) 

Here, the symmetric traceless product is defined by 

= Im< B j + A i B i) L_ 5 ij A h B k , (3.7) 

and the antisymmetric product is defined by 

A [i, B j] = l(A i B j - A j B l ). (3.8) 
For purposes of our calculation, it is also useful to define 

<o|0 o "( 3 p[ 8 i)|o> = £ (o\o»epf ] )|0) 

J=0,l,2 

= (Olx^-iSjV"^ Ptf ^ ( _|5 ) r (T r lr a x | 0) _ ( g g) 

It was shown in Ref. [26] that, by making use of the NRQCD equations of motion, one 
can express the LDME (0|e>f 3 ( 3 sj 1] )|0) in terms of the LDMEs (0|Of 1 ( 3 5[ 1] )|0) and 
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{0\O^ 2 { 3 S[ J )|0). (Equivalently, one can eliminate the LDME {0\O£ >3 ( 3 Sf l )\Q) by mak- 
ing use of a field redefinition [27].) Hence, we need not consider (0|O^ 3 ( 3 S'[ 1 ')|0) in our 
analysis. 

In the LDMEs (OlO^ 3 s[ 1] )|0> and (0|O^ 2 ( 3 5[ 1] )|0), one can replace \H+X){H+X\ in 
the projector T > h(m,o) with \H)(H\ (vacuum-saturation approximation), making an error of 
relative order v 4 . If one takes this approximation and evaluates the LDMEs in dimensional 
regularization in a potential model, then they are equal [28]. Since the static potential 
model is valid up to corrections of order v 2 , we have 

(0\O^( 3 S?)\0) = <O|0|f 2 ( 3 si 1] )|O> + Q(v 2 ). (3.10) 

Hence, up to corrections of relative order v 2 , only the sum of short-distance coefficients 
d4,i( 3 S^)(z) + d^ : 2( 3 S^)(z) appears in the fragmentation function. 

3.1 NRQCD factorization formulas for g — > J/ip through order v 4 

In summary, we have the following NRQCD factorization formulas for gluon fragmentation 
into J/ip through relative order v 4 . 
In relative order v° we have 

D [g -> J/J,] = d [g -> QQ(*sP)](0\Oi'*( 3 SP)\0). (3.11) 

The short-distance coefficient do[g — > QQ( 3 S^)] was calculated in Refs. [19, 20]. 
In relative order v 2 we have 

D 2 [g -+ JH>\ = d 2 [g QQ( 3 SP)](0\0 J 2 /Ip esn)\0). (3.12) 

The short-distance coefficient d 2 [g — > QQ( 3 S^)] was calculated in Ref. [21]. 
In relative order v 3 we have 

D 3 [g -> J/V] = d [g -> QQ( 1 4 8l )](0|O J/ ^( 1 4 8l )|0). (3.13) 

The short-distance coefficient do [g — > QQ^S^)] was calculated in Ref. [29] and differs 
from the short-distance coefficient do[g — > QQ( S )] in Eq. (3.13) only by a color factor, 
which we provide in Sec. 7. 

In relative order v 4 we have 

D 4 [g -> J/J,] = { d AA [g -> QQ( 3 S?)} + d A , 2 [g -> QQ( 3 S?)} } (OlO^sf^O) 
+ QQ( 3 P [8] )](0|O J/,/ '( 3 p[ 8 ])|0) 

+ -> 00( 3 ^P ] )](0|O^( 3 5p ] )|0). (3.14) 

The short-distance coefficient do[g — > QQ( 3 sf^)] was calculated at LO in a s in Refs. [21, 29] 
and at NLO in a s in Refs. [30, 31]. We verify the LO calculations in Refs. [21, 29] in the 
present paper, giving our result in Sec. 7. We compute d [g ->• QQ{ 3 P [S] )} in this paper, 
giving the result in Sec. 7. The short-distance coefficient do[g — > QQ( 3 P^)] was calculated 
in Refs. [29] and differs from the short-distance coefficient do[g —> QQ( 3 P^)] in Eq. (3.14) 
only by a color factor. The computation of the combination of short-distance coefficients 
^4,ib — > QQ{ 3 S^)] + d^2[g — > QQ{ 3 S^)] is the main goal of this paper. The result of 
that computation is given in Sec. 7. 
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4 Kinematics 



In the calculations to follow, in both full QCD and NRQCD, we employ the following 
kinematics. 

We take the Q and the Q to be free (on-shell) states with momenta 

P=lP + q, (4.1a) 
p=\P-q, (4.1b) 

respectively. The heavy quark has three-momentum q in the QQ rest frame, and, so, the 
invariant mass of the QQ state is 

P 2 = M 2 = 4E 2 , (4.2) 

where 

E = ^m 2 + q 2 . (4.3) 

We work in the frame in which the transverse momentum of the QQ pair vanishes. 
In this frame, the initial-state gluon, the final-state QQ pair and the final-state gluons, 
respectively, have the momenta 

k -(*,r ( , 4a) 



2k+ 
M 2 \ 

ki= (^ + '^T> fc ^)' ( 4 - 4c ) 

k 2 = [z 2 k + ,^-,k 2± ) , (4.4d) 



2z 2 k+ ' 

where we have introduced the longitudinal momentum fractions 



P+ 

z = — , (4.5a) 



Zl = ^, (4.5b) 



z 2 = ^. (4.5c) 

Because of the conservation of four- momentum, k = P + k\ + k 2 , the momenta k, k\ and 
k 2 depend implicitly on P, and, therefore, on q. We can make the dependence on q explicit 
by writing quantities in terms of dimensionless momenta 

- P 



It is also useful to express the Lorentz invariants in terms of the following dimensionless 
variables: 



ex = h ■ P, (4.7a) 
e 2 = k 2 ■ P, (4.7b) 
x = ki ■ k 2 = eie 2 (l - fci • k 2 ), (4-7c) 

where ki is the unit vector that is parallel to the three-vector ki in the QQ rest frame. 
The phase space in Eq. (2.7) can be expressed in terms of the dimensionless variables 

as 

4ttM 

d$ = -y^5(l-z), (4.8a) 
AnM d - 1 dzx d d ~ 2 k 1± 

4nM 2d - 3 d Zl dz 2 d d - 2 h± d d - 2 k 2± 

^ = 5AF e{2l) 6{Z2) *(1 - * " * " *) i^i^^)^^)^- ( 4 ' 8c ) 

We have not replaced the overall factor l/k + in Eq. (4.8) with l/(k + ^fP 2 ), because it 
ultimately will be cancelled by the factor k + in Cf rag in Eq. (2.5). The ranges of the 
variables z, z\ and z 2 are completely determined by the 5 and 9 functions. When we 
expand the fragmentation function in powers of q, it is convenient to make use of the phase 
space at LO in q, 

d$ n = d$ n \ q _ ¥0 , (4.9) 

where q — > means that, in the phase space in Eq. (4.8), we replace M with 2m. Then 
d& n can be expressed in terms of d<& n as follows: 

d$ = ~d$ , (4.10a) 
m 

/ F \ 3_2e 

d$x = I — J d$x, (4.10b) 
\m J 



d$ 2 = _ d$ 2 . (4.10c) 
\m J 

We express our results for the fragmentation contributions in terms of integrals over the 
phase spaces d<t n : 

D[9 -+ QQ] = D[g ^ QQ]\ mn ^ n . (4.11) 
The factors of E/m in Eq. (4.10) are then an additional source of relativistic corrections. 



5 Full-QCD calculations 

In this section, we compute the relevant fragmentation functions for free QQ states in full 
QCD. We have carried out the calculations by writing independent codes using reduce 
[32] and using the feyncalc package [33] in mathematica [34]. At each stage of the 
calculations we have checked that the independent codes give identical results. 
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The computations are carried out in d = 4— 2e dimensions with dimensional-regulariza- 
tion scale fi. We use the modified-minimal-subtraction (MS) scheme throughout. Then, 
in d = 4 — 2e dimensions, there is a factor [fj? exp(7 E )/(47r)] e that is associated with each 
factor of the strong coupling a s , where 7 E is the Euler-Mascheroni constant. 

In computing the QQ fragmentation functions, it is convenient to make use of projec- 
tion operators for the spin and color states of the QQ pair. The projection operators for a 
QQ pair in the color-singlet and color-octet configurations are 

Ai = -=1, (5.1a) 
Ag = V2T a , (5.1b) 

where 1 and T a are the identity matrix and the generator of the fundamental (triplet) 
representation of SU(3), a is the adjoint-representation color index (a = 1, . . . , — 1), and 
N c = 3. Spin-projection operators at LO in v were first given in Refs. [35-39]. Projectors 
accurate to all orders in v were given in Ref. [26]. For the spin-triplet state, the projection 
operator, correct to all orders in v, is 

P — i— *2E 

A(P,q,e* s ) = N(^-m) ^f-^—^ + m), (5.2) 

where eg is the spin polarization of the QQ pair, and N = [2y/2E(E + to)] -1 . Note that 
we use nonrelativistic normalization for the heavy-quark spinors. 

The use of the spin projection (5.2) in d dimensions requires some justification. It 
accounts for only the d — 1 vector polarization states, which, in the QQ rest frame, corre- 
spond to the d— 1 Pauli matrices Ui. In general, in d dimensions, one must consider states 
that correspond to products of the U{ that are linearly independent of the cxj [29]. These 
additional states vanish as e goes to zero. Hence, they can contribute only in conjunction 
with a pole in e. The poles in e in our calculation correspond to soft divergences. The 
divergent parts of soft interactions arise from the convection current on fermions lines, and, 
hence, do not change the fermion spin. Therefore, the additional states that correspond to 
products of the Oi never mix in our calculation with the vector states that correspond to 
the cij. Consequently, we need consider only the d— 1 vector states in our calculation. 2 

The spin-triplet, color-singlet part of an amplitude C is 

A4 = T!r[C(A<g>Ai)], (5.3) 

and the spin-triplet, color-octet part of an amplitude C is 

M a = Tr[C(A® Ag)], (5.4) 

where the traces are over the Dirac and color indices. The amplitude C includes the 
propagator of the initial gluon, as well as the associated polarization factor in Eq. (2.6). In 
our calculation, the amplitudes C, Ai, and Ai a are all expressed in terms of dimensionless 



2 Some elements of this argument were presented in Ref. [40] 



- 11 - 



variables in Eq. (4.6) or invariants that are formed from them, and so the dependence on 
q is explicit. 

The 5-wave part of M. (with color index suppressed in the color-octet case) can be 
written as an expansion in powers of v 2 = q 2 /m 2 : 



Ms = Mso + Ms2+M S 4 + 0(q b /m b ), (5.5) 



where 



Mso = (M),^ , (5.6a) 
2!(d-l) \dq a dqP J , 



MS2 = ^—r, I aP tt-^ , (5.6b) 



Mn- ^ r#* ( ^ \ (5 6c) 

54 4!(d - l)(d + 1) \dq<*dqPdqWq s J q ^ ' 1 j 



and 



jafl = _ g aP + p" P P (5.7 a ) 
ja/?7<5 _ ja/3pS _j_ ja-y j/35 _|_ ja8 j/3-y /g ^\ 

In order to project out the P-wave part of the amplitude .M, we multiply .M by the P- 
wave orbital-angular-momentum state — yJ~d^le* L ■ q and average over the direction of q. 
Here, ex, is the polarization vector for the orbital-angular-momentum state, and q = (0, q) 
in the rest frame of the pair. Then, the P-wave part of the amplitude is 

M P = Mpi + 0(q 3 /™ 3 )> (5.8) 

where 

We define squared amplitudes for the color-singlet and color-octet states as 

,4(*+i Z W) = C frag |^( 2s+1 /W)| 2 , (5.10a) 
A( 2s+1 lf) = C frag ^|^ a ( 2s+1 /f)| 2 , (5.10b) 

a 

where Cf rag is given in Eq. (2.5), and it is implicit that there are sums over the spin and 
orbital-angular-momentum polarizations of the QQ states and sums over the polarizations 
of the initial and final gluons. Note that 

£eT(A)ef (A) = ^ e L*W4(X) = I aP - (5.11) 

A A 

We denote the order-v fe contribution to A by Ak- 
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Figure 2. Feynman diagram for the fragmentation process g — > QQ( 3 S[ ). The dashed line 
represents the final-state cut. The momenta for Q and Q on the left side of the cut are p = \P + q 
and p = \P — q, respectively. The momenta on the right side of the cut are \P + q' and \P — q' , 
respectively. Here, \q'\ — \q\ in the rest frame of the QQ pair, but we distinguish the directions of 
q and q' in order to be able to project out orbital-angular-momentum states in the amplitude and 
its complex conjugate. 

5.1 D [g^ QQ( 3 sf)] 

The diagram for gluon fragmentation to a 3 Si color-octet QQ pair at order v° and a\ is 
shown in Fig. 2. A straightforward computation yields 

a s k + ( fi 2 



Carrying out the trivial integration over the phase space m Eqs. (4.8a) and (4.10a), 
we obtain 

D [g -+ QQ( 3 sf ] )\ = ^ (^ e7E ) <*(1 " *)■ ( 5 - 13 ) 
5.2 D 2 [g^QQ( 3 PW)) 

The diagram for gluon fragmentation to a 3 pl 8 l color-octet QQ pair at LO in a s and v is 
shown in Fig. 3. We obtain 



U ' 8N c (d-l)(d-2)m 8 \Att J ^ e?(l + 2ei) 2 ' [ ' 

where p n (z) are given by 



p (z) = (3 - 4e)(2 - 2e - 4z + 4z 2 ), (5.15a) 
Pl (z) = 2 [5 - lOe + 4e 2 - z(5 - 12e) + 2z 2 {l - 4e)] , (5.15b) 
p 2 ( z ) = 3 - 12e + 4e 2 + 2z(3 + 4e) - z 2 (5 + 4e), (5.15c) 
p 3 (z) = -2(1 - z) 2 . (5.15d) 
In some calculations in NRQCD, the P-wave orbital-angular-momentum state is normalized as — e* L ■ q. 
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Figure 3. One of the four Feynman diagrams for the fragmentation process g — > QQ( 3 P^). 
Three additional diagrams can be obtained by permuting the gluon-fcrmion vertices on the left and 
right sides of the cut. 



We carry out the integration over the phase space in Eqs. (4.8b) and (4.10b) by making 
use of the methods that are described in Appendix A. Then, we obtain 



D 2 [g^QQ( 3 pW)] 



8a 2 q 2 N 2 - 4 



(d - l)m 5 AN, 

-JL $(!-*) + /(*) 



(l-e)r(l + C ) 



E. 

47r' 



»7b 



Am 2 



p7e 



where the finite function f(z) is defined by 



1 - z) 1+2e 



+ 



4(1 



(l-z)-^-(l-z) 



-2e 



'13 - 7z) 



+(1 - zy 2e [10 + 3z- 5z 2 + 2e(2 - 2z + zf) - Ae 2 



+-(1 - z)- e [-28 + 152 + e(8 - llz) + Ae 2 z] 



Here, the distribution [g(z)] + is defined by 

/ 

J o 

In extracting the pole in Eq. (5.16), we have made use of the identity, 



dzh(z)[g(z)}+= / dz[h(z)-h(l)]g(z). 



1 



1 



(1-z) 



l+ne 



ne 



8(1 -z) + 



1 



(i - zy+ ne 



(5.16a) 



(5.16b) 



(5.17) 



(5.18) 



which applies when the domain of integration is < z < 1. The expression in Eq. (5.16b) 
gives the exact e dependence. We can expand the plus function [1/(1 — z) l+ne } + as 



.(1-*) 



l+ne 



E 

k=0 



(—ne) h 
kl 



log fc (l ~ z) 
1 - z 



(5.19) 
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Figure 4. One of the 36 Feynman diagrams for the fragmentation process g ->■ QQ( 3 sf ] ). Thirty- 
five additional diagrams can be obtained by permuting the gluon-fermion vertices on the left and 
right sides of the cut. 

In the analysis of f(z), we need to keep only terms through order e . Then, we can 
simplify f(z) as follows: 



/(*) 



1 



1-z 





"log(l - z) 




1-z 



1 



+ - [- 8 + 21 z - Wz 2 + 2(13 - 7z) log(l 



+ - [(23 - 16z)z + 5(8 - llz + 4z 2 ) log(l - z) - 3(13 - 7z) log 2 (l - z)] + 0(e 2 ). 
8 

(5.20) 

5.3 D 4 [g^QQ( 3 sP)} 

The diagrams for gluon fragmentation to a 3 Si color-singlet QQ pair at LO in a s are shown 
in Fig. 4. Through relative order v 4 , the relevant squared amplitudes are 



c^ s \M es?)\ 2 , 

2C iTag Re[M2( 3 S?)MU 3 S?)}, 

C frag (\M 2 ( 3 S?)\ 2 + 2Re[M 4 ( 3 sP)M* ( i s\ L] )} 



(5.21a) 
(5.21b) 
(5.21c) 



The order- v° contribution „4.o( 3 si^) ana - the order-v 2 contribution ^^S^) have been 
computed previously in Ref. [21]. Here, we wish to compute the order-u 4 contribution 

The integration of Ai^/S^) over the phase space contains soft divergences. These can 
arise when one or both of the final-state gluons become soft. We identify the divergent 
part of the integrand in Eq. (5.21c) that arises when both gluons become soft by making 
the substitutions k\ — > k\X and k 2 — > k 2 X, multiplying by A 4 , and taking the limit A — > 0. 
The result is 



5l2 



where 



C s (d-2) 

4 4 

e 1 e 2 
k + Tr 2 a 3 eZ d ~ 3 



[(d-2)e\el -2eie 2 x + 



N 2 



(d-l) 2 (d-2)m 8 AN 2 



An' 



=,7e 



3e 4 

SI 

m 4 



(5.22) 
(5.23) 
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We identify the divergent part of the integrand in Eq. (5.21c) that arises when only k\ 
(k 2 ) becomes soft by subtracting S\ 2 , making the substitution k\ — > k\\ (k 2 — > k 2 X), 
multiplying by A 2 , and taking the limit A — > 0. The result is 



Si 



e\e\(l + 2e 2 Y 



y^x n h n (e 1 ,e 2 ,z 1 ,z 2 ), 



n=0 
2 



S2 efef(l + 2 ei ) 2 



(5.24a) 
(5.24b) 



ra=0 



where h n (ei, e 2 , z\, z 2 ) are given by 



ho(ei, e 2 , z 1 ,z 2 ) = 2e 2 (l + 2e 2 ){ - 2(3 - 2e)(l - e)e 2 e 2 + zieie 2 + z 2 e?e 2 

x [1 - 2(1 - 2e)e 2 ] - z 2 e^(l + ee 2 ) + ziz 2 eie 2 [l - 2(1 - e)e 2 ] 
+z 2 2 e 2 [-l + (1 - 3e)e 2 + 2(1 - 2e)e 2 , ] } + 4(3 - 2e)(l - e)e 2 e|, 

(5.25a) 

hi(ei,e 2 , zi,z 2 ) = 2e 2 (l + 2e 2 ){zie 2 (l + 2ee 2 ) + z 2 ei[l + (3 - 2e)e 2 ] - Z\Z 2 {1 + 2ee 2 
+2eel) - z|ei[3 - 2e + 2(1 - e)e 2 ]} + 2eie 2 .[5 - 4e + 6(1 - e)e 2 ], 

(5.25b) 

h 2 (e 1 ,e 2 ,z 1 ,z 2 ) = -(1 + 2e 2 ){e 2 (4 - 3e + 2ee 2 ) + 2z 2 e 2 (l - 2e - 2ee 2 ) 

-2z 2 (l + e 2 )(l - e - ee 2 )} - (2 - e)e 2 . (5.25c) 

We carry out the integrations of S± 2 , S± and S 2 over the phase space d<& 2 in Eq. (4.10c) 
[see also Eq. (4.8c)] by making use of the methods that are described in Appendix A. Then, 
we obtain 

1 - z 1+2t 



I[S 



12 



1 

8^ 



8(1 -z) 



1 

2em 



1 

{l^z) 



l+4e 



+ 



2ei R (l - z) l +^ 



V37rm 2 y 16iV 2 (d-l)m 3 V4vr / d-1 



I 4m 2 



fj^^V r 2 (i + e)r 2 (i-2 £ ) 



T(l-4e) 



(1 - e)(6 - 2e - e 2 - 2e 3 



(5.26) 



where we have used the identity in Eq. (5.18). If we expand Eq. (5.26) in powers of e, then 
we find that 



I[S 



12 



1 



5(1 - z) 



1 



+5(1 - z) 
+2 



2e m 
1 - 3vr 2 2 



5(1 -z) 



- - log — I + 

3 6 2m 



1 



1 - js 



48 



- log — + log 2 — 
3 2m 2m 



+ 



1 - z 



1 



2 u 
- - 2 log — 

3 & 2m 



log(l - z) 
1 - z 



21 : -iog^ r V4 lo ^-2iog(i-.) 



8a, \ 2 A^ 2 -4 



2m / 1 — z 
6q 4 



V3vrm 2 y 16iV c 2 (d - l)m 3 \An^ J d-1 



+ 0(e) 



(5.27) 
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We also find that 

J[Si] = I[S 2 ] 



n 

2 M2 



\4m 2 



2e Z- 2 + 2e 



/ 8a, y iV 2 -4 (n 



7e 



7ra s 



(l-z)~ 4 T 2 (l + e 
48(1 - e) 

+ 0(e), 



V3vrm 2 y 16iV| V 47r / (d - l)m 3 (f - 1 
where 

r = 3 j 12(13 - 7z)z 2 Li 2 (l - z) + 3(7z - 13)z 2 log 2 (l - z) 

-2z [z 3 - (43 + 7tt 2 )z 2 + 6(2z + 9)z 2 log z + (48 + 13vr 2 )z - 42] 
+3(8z 4 - 7z 3 + 34z 2 - 44z + 28) log(l -z)\, 
n = 18z 2 [z(21 - lQz) + 2(13 - 7z) log(l - z)\ . 
Here, Li 2 (;z) is the Spence function, which is defined by 

log(l - t) 



(5.28) 



Li 2 (x) 



dt 



OO U 



£ 

fc=i 



A; 2 



The fragmentation-function contribution of ^4( 3 S'{ 1 ') is then 



DAa -»• Q0( 3 M 1] )] = Ai[s ->• QQ( 3 M 1] )] finitc + i[S 12 ] + i[Si] + i[s a 



where 



Ab^QQ( 3 s! 1] )] finitc 



(5.29a) 
(5.29b) 

(5.30) 

(5.31) 
(5.32) 



Since D[g — > QQ( 3 5[ 1 ')] finite contains no soft divergences, we compute it in d = 4 dimen- 
sions, using numerical integration over the phase space. This computation is described in 
Sec. 8. 

As was mentioned previously, D±[g — > QQ( 3 S^)] &mtc also contains contributions that 
arise from the difference between the fully relativistic phase space (i$ 2 and the order-v 
phase space d<3? 2 [Eq. (4.10c)]. Hence, we write 

D±[g -»• QQ( 3 S?)\ = D 4 [g QQ( 3 sP)f n[tc + I[S 12 ] + J [Si] + I[S 2 ], (5.33) 



where 



5q 2 



3 C [1]^ 



D^g -> gO( 3 M 1] )] finite = Atb -> QQ( 3 M 1] )] finite + f^^b -> W( 3 s| 

2m z 



where 



D [g^QQ( 3 S?)}= I d4> 2 Ao( 6 S l i l ), 



D 2 [g^QQ( 3 sf ] )] = I d^AtCS^), 



C3 [l]> 



(5.34) 

(5.35a) 
(5.35b) 
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and we have used the expansion 



Do[g — > QQ( 3 S^)] and D2[g — > QQ^S^ 1 )] are finite, and they have been evaluated in d = 4 
dimensions in Ref. [21]. We have checked those calculations, computing the phase-space 
integrations in D [g -> QQ( 3 s[ 1] )] and D 2 [g -> QQ( 3 flf ] )] numerically. 



(5.36) 



6 NRQCD LDMEs 



In this section we compute the NRQCD LDMEs for free QQ{ 3 S^) states that are relevant 
through relative order v 4 . These computations are carried out in each case at the leading 
nontrivial order in a s in d = 4 — 2e dimensions with dimensional-regularization scale fi. We 
remind the reader that, because we use the MS scheme in computing the QCD corrections 
to NRQCD LDMEs, there is a factor [/x 2 exp(7 E ) / (47r)] e that is associated with each factor 
of the strong coupling a s in d dimensions. 

6.1 Order a° s 

The matrix elements of the QQ NRQCD operators at order oP s are normalized as 

(6.1a) 
(6.1b) 
(6.1c) 
(6.1d) 
(6.1e) 
(6.1f) 

where a sum over the final-state polarizations is implied. The superscript (k) indicates the 
order in a s . 

6.2 Order a s 

6.2.1 3 P® -> 3 s[ 1] 

In order a s , the 3 pl 8 l QQ operator can couple to the state through the diagrams that 
are shown in Fig. 5. 4 

The set of diagrams in Fig. 5 is gauge invariant. We find it convenient to work in 
the QQ center-of-momentum frame and to compute the diagrams in the Coulomb gauge. 





|o}(°) 


= {Nl 


-1), 








|o}(°) 


= 2{d- 


- i)n c , 






(o|0? O(3 ^ 81) ( 3 sf ] ; 


|o)(°) 


= (d~ 


" 


1), 




(0 | O w^ b] )(3 5 n : 


|o)(°) 


= 9 2 (0 




] ) ( 3 5 [n 


)|o><°> 


(o|of (35[n]) ( 3 M" ] ; 


|o)(°) 


= 9 4 (0 




] ) { 3 s [n 


)|o><°> 


(o\of {3pls]) ep®: 


|o>(°> 


= q\d 


- l)(^c 2 







4 We suppress Wilson lines in diagrams involving NRQCD production operators [3, 4] because the dia- 
grams involving interactions with Wilson lines vanish for the orders in a s and the operator and final-state 
quantum numbers that we consider. 
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Q(E,q) Q{E, q') 

Q(-E,q) Q{-E,q') 



Figure 5. One of the four Feynman diagrams for the computation of the LDME 

(0|C^ (5(3 ^ 11) [OQ( 3 P[ 8 ])]|0)( 1 ). The solid circles represent the QQ operators in the LDME. As 
in the full-QCD calculation, we take the free Q and Q states to have momenta \P + q and \P — q 
on the left side of the cut and \P + q' and \P — q' on the right side of the cut, where, |q'| = \q\ in 
the rest frame of the QQ pair, but we distinguish the directions of q and q' in order to be able to 
project out orbital-angular-momentum states in the amplitude and its complex conjugate. Three 
additional diagrams can be obtained by permuting the gluon-fermion and operator vertices on the 
left and right sides of the cut. 



Then, the real gluons in the final state must be transverse. Through relative order v 2 , a 
straightforward computation in dimensional regularization gives 

( 0| W^)(»pM)|o>a) = (mU+mP+mV+M®) , (6.2) 



where 



d^k (k + 2q) l (k + 2q) i {6 i i - + 2q'y(k + 2q') 1 



2| fc| (2vr) rf — 1 {k 2 + 2k- q + 2m\k\)(k 2 + 2k- q' + 2m|fe|) ' 

f d d ^k (-k + 2q) l (-k + 2q) i {5 i i -k i ti){k + 2q'y(k + 2q') 1 
X J 2|fc|(2vr) a! - 1 (k 2 - 2k ■ q + 2m\k\){k 2 + 2k ■ q' + 2m\k\) 

d d - x k {-k + 2q)\-k + 2q) i {5 ij - fc l ^')(-fc + 2q'y(-k + 2q') 1 



(6.3a) 



(6.3b) 



2|fc|(2vr) d - 1 (k 2 - 2k ■ q + 2m|fc|)(fc 2 - 2k ■ q' + 2m\k\) 

d d - x k (k + 2q) l (k + 2q) i {8 lj - + 2q'y{-k + 2q') 1 

2|fc|(27r) ci - 1 (k 2 + 2k ■ q + 2m\k\){k 2 - 2k ■ q> + 2m\k\) ' 



(6.3c) 



(6.3d) 
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In Eq. (6.3), £ and are the Pauli spinors for the free Q and Q states, respectively. The 
subscript 3 S , [^ in Eq. (6.2) indicates that the bispinors ^r] and 77^ are in color-singlet, 
spin-triplet states and that we project onto S- wave states by averaging over the directions 
of q and q' . A sum over polarizations of the spin-triplet QQ pair is implicit. 

We expand the integrands in Eq. (6.3) in powers of 1/m. In dimensional regulariza- 
tion, only the leading power contributes because the expressions for higher powers in 1/m 
produce power-divergent, homogeneous integrals. 5 The result is 



35W 



(0 | O f ( 3 4 1] ) ( 3 P[ 8 ])|0)(1) = a™. d_^2 /|-i ( cWVWgV ^ 

|fe| 3 (27r) d - 1 ' { } 

The remaining integration over the spatial components of k is straightforward and 
yields 

<o| r si '\^)|o>« - - ^ (eWAiW*)^ . 

(6.5) 

Here, we have separated the ultraviolet (UV) and infrared (IR) divergent contributions of 
the scaleless integral. The quantity c(e) is given by 

m (^r(i-6)r(i) 

(1-36) (l-2e)r(i-e) 



Note that 



Now, 



c(0) = 1. (6.7) 



£<!l = J- + *>Zi, (6.8a) 
euv euv e 

em em e 

where we have dropped the subscripts "UV" and "IR" in the second terms of the above 
equations because those terms are finite. Hence, we have 



( ) ( 1 1 ^ - 1 1 (6 9) 

2euv 2eiR/ 2euv 2eiR 



Therefore, 



(o| 0o ^'Y P , 81)|0)(1) = ^ - 4) (^A,V^f): 



(6.10) 



J This approach was first used in Appendix B of Ref. [2]. It has been discussed subsequently in Refs. [29, 
41, 42]. In Ref. [43], it was pointed out that this approach allocates contributions that are infrared finite 
to the short-distance coefficients, rather than to the LDMEs, and is, therefore, the NRQCD analogue of 
the standard methods for computing dimensionally regulated short-distance coefficients for hard-scattering 
processes in collinear factorization in QCD. 
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We renormalize (0|O^ (35 ' ]) ( 3 P [8] )|0) (1) in the MS scheme. In the MS scheme, one 
renders the loop corrections to operators finite by subtracting pure poles in e\jy, setting 
d = 4 in the loop-correction factors that multiply the poles, e-dependent factors that are 
not associated with the loop correction, such as those that arise from the projections of 
operators onto physical states, are evaluated in d dimensions. Hence, we find that the MS 
counterterm is 



35W 



Z (i) ( 3 p[ 8] ^ 3^1]) _ x ] ^WaSvW^)^ . (6.11) 



Then, we have 



(o\of {3sl i ] \ 3 pW)\o) ( A = (o\of (is i ] \ 3 pW)\o)W + 5( |o Q<5(3sill) ( 3 p[ 8 ])|o)( 1 ) 



MS 

8a, / -1 \ N? - 1 



CVVA^WA) m- (6-12) 



3im 2 VW 4iV| V s * * '''' * Vs S K 

Now we extract the 5-wave part by averaging over the angles of q and q': 

f ! s! 1 h,,j«K,.,m 8« s / -1 \ iVg - 1 1 ,„ 1OT oo^W) 
37rm 2 V 2 eiR/ 47V C 2 d 

(6.13) 



6.2.2 3 5p ] -> 3 Pl 8 ] 

In order a s , the 3 S® QQ operator can couple to the QQ( 3 P^) state. By carrying out a 
calculation that is very similar to the one in the preceding section, we find that 

ZiDfsf] 3p[8]) _ 1 = p, (j±\ (6 14) 

v 1 ; 3vrm 2 V2euv / 4iV c v ' 

and that 

<0|^>flf )|0>« = £± (^) ^(0|^^ [8l ) ( 3p[«I)|0)(°). (6.15) 

The equivalent result for NRQCD decay LDMEs was obtained in Ref. [40]. 
6.3 Order a 2 s 

In order a 2 s , the 3 5® operator can couple to the QQ( 3 S^) state through the diagrams 
that are shown in Fig. 6. As in the order-a s case, we expand the integrands for these 
diagrams in powers of 1/m and, again, only the leading power contributes in dimensional 
regularization. Now we combine contributions that differ in the order of the gluon vertices 
on a quark or an antiquark line by making use of the fact that the color factor for the 
color-singlet contribution is symmetric under the interchange of the gluon color indices 
and by making use of the identity 

111111 , . 

+ — ^ = mu\- ( 6 - 16 ) 



\k\ \k\ + \£\ \£\ \k\ + \£\ \k\ \£\ 
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Q(E,q) 



Q(E, q') 




Figure 6. One of the 36 Feynman diagrams for computing the LDME 

(0|O^ (5(3 ^ 11) [0Q( 3 S'[ 8l )]|0) (2) . The solid circles represent the QQ operators in the LDME. Thirty- 
five additional diagrams can be obtained by permuting the gluon-fermion and operator vertices on 
the left and right sides of the cut. 



After we combine the contributions of all of the diagrams in this way, the loop integrations 
decouple, and we have 



<o|of (35i >sl 81 )|o> {2) 



(4vra s ) 2 (tf 



4tt 



2e (iV c 2 -4)(^-l) 



Jd-lk $U - f d d ~ l l 5 rs 



(2Tr) d - 1 \k\ 3 
8a s ( 1 1 



ft 



(27T) 



d-1 



3irm 2 \2euv 2erR 



(JV*-l)(JV*-4) 



16iV 3 



*1 



(6.17) 



where we have made use of Eq. (6.9). 

We carry out the renormalization of {0\O^ 3Sl ) ( 3 Sf ] )\Q)W in the MS scheme. First, 
we add the contribution of the one-loop diagrams involving the counterterm ^( 1 )( 3 p[ 8 l — y 
3 S[ 1] ) - 1, which are shown in Fig. 7. The contribution of these diagrams is 



5 1 (0\of {3sl ' ] \ 3 sfW 2) 



where (0\Oq Q( - 3s ^\ 3 P^)\0)^ is given in Eq. (6.10) and Z^{ 3 sf ] -»■ 3 p[ 8 l) - 1 is given 
in Eq. (6.14). Thus, 

(JV*-l)(jV*-4) 



Z (1)(3 5 J8] ^3p[8] ) _ 1 



| O W( 3 4 1I ) ( 3 p [8] ) | 0) (l) ) (61g) 



s 1 (o\of i3sl ' ] \ s sfW 2) 
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Q(E,q) 



Z( 3 S IS] ^3p[8] ) _ 1 



Q(~E,q) 



Q(E,q') 




Figure 7. One of the four Feynman diagrams for computing the one-loop correction to the coun- 
terterm Z( 3 sf ] 3 pl 8 l)-l. The symbols (g) represent the counterterm Z( 3 sf ] 3 pl 8 l)-l. Three 
additional diagrams can be obtained by permuting the gluon-fermion and counterterm vertices on 
the left and right sides of the cut. 

Adding this counterterm contribution to Eq. (6.17), we obtain 



|0 QQ( 3 ^) ( 3 5 J8] )|0)(2) 



1 / 8a s V 
8 V3vrm 2 J 



■1 



+ 5 1 (0\of {3S ^\ 3 sfW 2) 

2 lV/\r2 



+ 



1_\ (N C 1)(N C 4) (fff^ff^ ((k , ()) 



, uv c IR/ 16tf| v /3sr 

We see that the counterterm contribution removes the cross term between the pole in euv 
and the pole in eiR. Hence, the remaining overall UV divergence is not coupled to infrared 
contributions, as is required by the consistency of the renormalization program. There are 
no single poles in euv- The remaining double pole in euv is removed by adding the MS 
counterterm 



5 2 (o\of {3sl ' ]) esfW 2) 



\3-7rm 2 / 8e uv 



16 N? 



£q H q' r <r n rivU i q r cr n ( 



3 S [1] 
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Adding this counterterm to Eq. (6.20), we obtain 

f3?['h o rci (2) nniS^U 

'MS 



|o QQ(^) ( 3 5 j8] )|0) g) = (0 | O W^) ( 3 5 M ) | 0) (2) +5l(0 |O W ^ 1, )(3 5 p] ) | 0) (2) 



+5 2 (0|O Q<5(3silI) ( 3 M 8] )|0) (2) 
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1 (iV2-i )( Ar2_4) 



37rm 2 / 8ef R 
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(6.22) 
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We evaluate the S- wave part by carrying out the average over the directions of q and q' , 
with the result 

«<3<»si"W|8h lm e> . 1 ( 8a, V (A? - 1)(JV C 2 - 4) 



(0| O ««^I")(3 s | 8 ] )|0) g) = 1 C 

° e IR V 



, |; v 37rm 2 J 16JV| 
x -i— (0| Ojf (35 ' 1] } ( 3 ,SJ 1] ) 1 0) (0) . (6.23) 

6.4 Renormalization-group evolution 

In the renormalized LDMEs, we can identify the dimensional-regularization scale fi with 
the NRQCD factorization scale fj,\. We now work out the renormalization-group evolution 

with respect to /x A of the LDMEs (0\o9 Q{3 ^ ( 3 P [8] )|0)^, (0\o9 Q{3p[S] \ 3 sf ,] )\0) { ^ and 
(0\O® Q{3Sll \ 3 sf ] )\0)® First, we take d/dlogfi A of Eq. (6.13): 



'MS 



dlo g/ u A N 1 u v ' MS dlogn A da s x 1 u v 'ms 

_ _8o^ N^l J_ QQ(3 S W) 3 [1] (o) 

" 3W 4iV c 2 d-l^i ( 61 j|U) ' 

(6.24) 

where we have used 

J a ' s = -2ea s + 0{a 2 s ). (6.25) 
dlog^A 

The result in Eq. (6.24) agrees with the corresponding results of Refs. [23, 44]. 
Similarly, taking d/dlogfi A of Eq. (6.15), we obtain 

^(0|0?«< >PM) ( 3 S! 8, )|0>S = 3^4^<°l°»' 5<,J " ! ' ) ( 3i " 81 »l > (0) ' (6 - 26 » 
Taking d/d\og[i\ of Eq. (6.23), we also obtain 

rfl^7^ ( ' ° C l}l ) SS _ lwJ V2l^y 16iV| 

x^l y (0|O 4 Q f (35 ' 1,) ( 3 5[ 1] )|0)(°). (6.27) 

Substituting Eq. (6.13) into the right side of Eq. (6.27) we find that 

" -(0|O?« N'VsfbKe = J^%i(0|O?« <3S !"»(»pM)|0,g). (6.28) 



dlog^A MS 3vrm 2 4iV c N 1 v n MS' 

Eqs. (6.26) and (6.28) agree with the result in Eq. (B19b) of Ref. [2] at the leading nontrivial 
order in v and with the corresponding result in Ref [44] , but disagree with the corresponding 
result in Ref. [23]. 
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7 NRQCD matching and short-distance coefficients 

In the discussion to follow, and in remainder of this paper, the short-distance coefficients 
d n that appear are always in the MS scheme. For brevity, we do not indicate the scheme 
explicitly. 

In relative order v 4 , the NRQCD factorization equation that relates full QCD and 
NRQCD for the fragmentation of a gluon into a 3 S^ QQ pair is 

D A [g -»■ QQ( 3 S?)p = {d^[g -+ QQ( 3 S?)}^ + d A , 2 [g -> QQfS?)}®} 

x(0|of (35[1]) ( 3 5l 1] )|0)(°) 
+ d [g -> QQ( 3 P[ 8 1)]( 2 )(0|0 QO(35[1I) ( 3 P[ 8 ])|0)« 
+ d [g -> QQ( 3 M 8l )] (1) (0|O Q<3(35il]) ( 3 5l 8] )|0)( 2 ), (7.1) 



where we have shown the contributions at the leading nontrivial order in a s , namely, a 3 , and 
the LDME (0\O^ Q{3S ' ] \ 3 sP)\0)^ means either (0\off ( 3 sf ] )\0)^ or (0|C>jf 
( 3 S , ( 1 ^)|0)(°), since they are equal at the present order of interest in v, as we have explained 
in Sec. 3. We wish to use this matching equation to determine the sum of short-distance 
coefficients d 4 ,i[9 ~> QQ^S^)]^ + d^ 2 [g -> QQ( 3 S^)]( 3 \ We have already computed the 
QQ LDMEs on the right side of Eq. (7.1). In order to fix the short-distance coefficients 
in the second and third terms on the right side of Eq. (7.1), we make use of a matching 
equation at relative order v 2 , which is 

D 2 [g -+ QQ( 3 P [8] )] (2) = d [g -+ QQ( 3 pW)fH0\of {3plS] \ 3 pW)\0)^ 

+ d [g -> QQ( 3 sf ] )]^(0\O^ 3pl8] \ 3 sf ] )\0)^, (7.2) 

and a matching equation at relative order v°, which is 

D [g -»• QQ( 3 sf ] )}^ = do [g -»■ QO( 3 M 81 )] (1) (0|0 Q<5(35 ' 8,) ( 3 S[ 8] )|0) (0) . (7.3) 

First, we solve Eq. (7.3) for d [g -> QQ^Sj 81 )] W, making use of Eq. (5.13) for D [g -> 
<2Q( 3 S{ 8] )]« and Eq. (6.1c) for (0|oJ 0(3 ^ 81) ( 3 5j 8] )|0)(°). The result is 

4b - <J0PfiP)l« - ( „_ 1)( %_ 1)ro , " •)■ C"> 

The short-distance coefficient in Eq. (7.4) agrees with that in Refs. [30, 31]. 

Next, we determine the short-distance coefficient do[g —> QQ( 3 P^)]^ by making 
use of the order-?; 2 matching equation (7.2). We substitute D 2 [g — >■ QQ( 3 P^)]^ in 
Eq. (5.16), d [g QQ( 3 sf)]W in Eq. (7.4), (0|Oj 0{3p[81) ( 3 p[ 8 l)|0)(°) in Eq. (6. If), and 
(O|0j O(3p[8 VM 8l )l°> (1) inE q- ( 6 - 15 ) intoEq. (7.2). Then, solving for d [g -»• QQ( 3 p[ 8 l)]( 2 ), 
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we obtain 

d [g^QQ( 3 P [8 ¥ 2) 



8ai 



3(d- l)(iV 2 - l)m 5 



AN, 



MA^E 

4ir 



5(1 -z) 
2em 



(1 - e)r(l + e) / /ij t 
\4m 2 



:,-/E 



+ 



(i- e) r(i + t) , 



V 4m 2 



(7.5) 



where /(z) is given in Eq. (5.16b). In Ref. [29], the color-singlet short-distance coefficients 
do[g — > QQ( 3 Pj were computed. Summing the results for do[g — >■ QQ( 3 Pj )]^ in 
Ref. [29] over J = 0, 1 and 2 and multiplying by [(iV 2 - l)/(4iV 2 )]- 1 [(iV 2 - 4)/(4iV c )] in 
order to obtain the corresponding short-distance coefficient for the color-octet channel, we 
find agreement with our result in Eq. (7.5). 

The expression in Eq. (7.5) gives the exact e dependence. Expanding this expression 
to order e , using the expression for f(z) in Eq. (5.20), we find that 
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(7.6) 



Finally, we determine the sum of short-distance coefficients d^\[g — > QQ{ 3 S^)]^ + 
d±, 2 [g -»• <2Q( 3 s[ 1] )] (3) by making use of the matching equation at relative order (7.1). 
We substitute D A [g -> QQ( 3 ,s| 1] )]( 3 ) in Eq. (5.33), do [5 -> QQ( 3 P [8] )] (2) in Eq. (7.5), d [ff -»■ 

gg( 3 M 81 )] (1) inE q- (7.4), (o\o® Qes i ]) ( 3 s?)\o}W m Eq. (e.ie), (o|0? Q ^V-P pl )|o) (1) 

in Eq. (6.13), and (Q\O% Q ? S ™\ 3 S^ ] )\0)W in Eq. (6.23) into Eq. (7.1). Solving for d 4>1 \g -> 
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QQ( 3 S?)}® + d Ay2 [g -> QQ( 3 5l 1] )]( 3 ), we obtain 
dt,x[g -> QQ( 3 M 1] )] (3) + ^4,2 [3 gQ( 3 M 1] )] {3) 
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where we have used the identity Li2(l 
Qg( 3 si 1] )] finite is defined by 



ir 2 /G - Li 2 {z) - logz log(l 



^b^Q0( 3 M 1] )] finitc 



D 4 [g^QQ( 3 M 1] )] finitc 
{o|cf 0(3 ^ 11) ( 3 s , l 1] )|o)(°) 



3 clilNlfinitc 



D,[g^QQ( 3 S[ Li )] 
2{d-l)N c q± 



(7.7) 
<k\g -> 

(7.8) 



where D A [g -> QQ( 3 5f ] )] finitc and (Olcf^'V^ 11 )! ^ are given in Eqs. (5.34) and 
(6.1e), respectively. The result in Eq. (7.7) is new. As expected, both the double and 
single poles in erR have cancelled in Eq. (7.7). These cancellations rely nontrivially on the 
correctness of the infrared subtractions and on the NRQCD operator renormalizations in 
our calculation. 

Finally, we note that the color-singlet short-distance coefficient do[g — > QQi 1 S^)}^ 
was computed in Ref. [29]. We can obtain the corresponding color-octet short-distance 
coefficient by multiplying by the ratio of the color-octet and color-singlet color factors, 
namely, [(iV c 2 - l)/(4iY 2 )]^ =3 [(iV 2 - 4)/(42V c )] . Then, we find that 



8m 3 4N C 



[3* - 2z 2 + 2(1 - z) log(l - z)] 



(7.9) 



8 Numerical results 



In this section we describe the numerical results that derive from our calculations. 

We have evaluated d^[g — > QQ( 3 S^)] 5mtc (z) in Eq. (7.8) by carrying out the integra- 
tions over the phase space ff<3?2 in Eqs. (5.32), (5.35a), and (5.35b) numerically at 100 points 
in each of the ranges z = to z = 10~ 2 , z = 10~ 2 to z = 1 — 10~ 2 , and z = 1 — 10~ 2 to z = 1. 
We have then used the parametrization in Eq. (B.2) to obtain a best fit to the numerical 
results, which leads to the parameters that are given in Table 3. Details of this procedure 
are given in Appendix B. We have also applied this procedure to do[g — > QQ( 3 S^)]( 3 \z) 
and d 2 [g —> QQ( 3 S^)]( 3 \z). In Fig. 8, we plot the results of the numerical calculations of 
the short- distance coefficients. 
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Figure 8. The color-singlet short-distance coefficients d^[g QQ( 3 S[ )]^ 3 \z), d%\g — > 
QQ( 3 S [ ^ ] )]^(z) and d t [g -> QQ( 3 S'[ 11 )] finito (z), which are defined in Eqs. (3.11), (3.12) and (7.7), 
respectively, as functions of z. The scaling factors are (Af , Af 2 , A/4) = (10~ 3 x a 3 /m 3 ,10~ 2 x 
af/m 5 , 10~ 2 x a 3 /m 7 ). 

We would like to obtain estimates of the relative sizes of the contributions of the 
various fragmentation functions to the cross section for J ftp production in hadron-hadron 
collisions. The contribution of a fragmentation process to the J/i/j production cross section, 
differential in the J ftp transverse momentum p T , is 



where ^~{Pt) ls the cross section to produce a gluon with transverse momentum p T . Let 
us assume that (da g /dp T )(p T ) oc 1/p^, where k is a fixed power. Then, 




(8.1) 



dp T 



oc I K {D) 



(8.2) 



where 




(8.3) 
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Table 1. Results of computing I K (d) for the order-?; , order- 1> 2 , and order-u 4 fragmentation func- 
tions for the fragmentation process g — > QQ^S 1 ^). The first two rows are the normalization inte- 
gral I Q (d) for the short-distance coefficients d [g -> QQ^sj 11 )]^ [Eq. (3.11)], d 2 [g -» Q<9( 3 -S'i 11 )] (3) 
[Eq. (3.12)] and d iA [g QQ( 3 S l ^ ] )]^ + d^\g -> QQ( 3 S{ ] )]W [Eq. (7.7)]. The third and fourth 
rows are I (d) times the ratio of LDMEs R n = (0\O n ( 3 S l ^ ] )\0) / (0\O ( 3 S l ^ ] )\0) , as is described in 
the text. The fifth and sixth rows are the integral 15.2(d)- The seventh and eighth rows are 15.2(d) 
times R n . In the case of dt,i[<? -> QQ( 3 S^)}^ + d^lg -> QQC 3 ^)]^, we have set fi\ = m and 
2m. The quantities in the table are divided by the factor AT = 10 -4 x a 3 /m 3 . 



Ud^Xdnig^QQfS^)]^ 


do 


d 2 


^4,1 + C?4,2 


Io(dn)\fj,A=m 


8.29 


20.3/m 2 


-37.4 /m 4 


Io(d n )\u A =2m 


8.29 


20.3/m 2 


-6.54/m 4 


Io(dn)\n A =m X Rn 


8.29 


3.95 


-1.41 


Io(d n )\ u A =2rn X Rn 


8.29 


3.95 


-0.247 


h.2(dn)\fM A =rn 


0.743 


11.2/m 2 


35.7 /m 4 


h.2(d n )\nA=2m 


0.743 


11.2/m 2 


84.5 /m 4 


h.2(d n )\ij, A =m x Rn 


0.743 


2.18 


1.35 


h>.2(d n )\i±A=2m, X Rn 


0.743 


2.18 


3.19 



Hence, we can obtain a rough estimate of the relative contribution of a fragmentation 
process to the cross section by computing I K (D). & 

We can estimate k by taking advantage of the fact that J/ip production in the 3 5^ 8 ' 
channel is dominated at LO in a s at large p T by gluon fragmentation into a J/ip with 
longitudinal- momentum fraction z = 1 (Ref. [45]). Consequently, the p T dependence in 
this channel at high p T is approximately that of da g /dp T . Furthermore, since the NLO 
k factor for this channel is essentially independent of p T and amounts to a correction of 
only about 14% at the Tevatron [46], we do not expect the NLO corrections to change the 
p T dependence of this channel significantly. Therefore, we estimate k by making use of the 
result for the 3 S*[ 8 ' contribution to do jdp T at NLO in a s that appears in Fig. 1(c) of Ref. [7]. 
Specifically, we compute (d/dlogp T )log[da/dp T ] and find that k « 5.2 at p T = 20 GeV. 

We have computed I K (d) for k = and 5.2, taking d to be the order- v°, order-u 2 and 
order-w 4 fragmentation functions for the fragmentation process g — > QQ( 3 SP). The results 
are shown in Table 1. In the first and second rows of Table 1, we give the normalization 
integrals Io(d) for the short-distance coefficients. In the third and fourth rows, we give 
Io(d) multiplied by the relative value of corresponding LDME. In the fifth and sixth rows 
we give 15.2(d). In the seventh and eighth rows we give 15.2(d) multiplied by the relative 
value of corresponding LDME. We take the relative values of the LDMEs to be given by 

6 The precise calculation of a fragmentation contribution to a cross section would require that one com- 
pute the convolution of the two-to-two partonic cross sections that produce a final-state gluon with the 
fragmentation function and with the appropriate parton distributions. Such a calculation is beyond the 
scope of the present paper. 
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the generalized Gremm-Kapustin relation [28]: 

(OIO^C^JIO) : (0\O^( 3 S?)\0) : <O|0 4 J/ ^( 3 s[ 1] )|O> = 1 : m 2 (v 2 ) : m A (v 2 ) 2 . (8.4) 
We use the value of m 2 (v 2 } from Table I of Ref. [47]: 

m 2 (v 2 ) = 0.437 GeV 2 at m = 1.5 GeV. (8.5) 

Examining the relative values in the fifth or sixth rows of Table 1, we see that, if we 
exclude the factors from the NRQCD LDMEs the order-?; 4 contribution enhanced consid- 
erably relative to the order-?; and order-f 2 contributions. We attribute this enhancement 
to the strong peaking of 04,1b — > QQ( 3 S^)]^ + 04,2 b — > QQ( 3 S^)]^ near z = 1, which 
arises primarily from the 5 functions and + distributions that it contains. We expect 
such strong peaking to occur whenever the full-QCD process contains soft divergences that 
correspond to divergences in the QQ LDMEs. 

From the seventh and eighth rows of Table 1, we see that the order-?; fragmentation 
contribution to the cross section is enhanced by approximately a factor of 2 for fi\ = m 
and by about a factor of 4 for fj,\ = 2m, relative to the order-?; contribution. However, 
the order-?; contribution to the J ftp cross section at p T = 10 GeV lies about a factor 
of 30 below the measured cross section. Thus, while the enhancement of the order-?; 
fragmentation contribution to the cross section is substantial, it is not sufficient to make the 
order-?; 4 fragmentation contribution to the cross section important for the phenomenology 
of J/ip production at the current level of precision. 

We can also compare the relative contributions to da^J^ of each of the three QQ 
channels in that appear in D[g — > J/ip] m order v 4 [Eq. (3.14)]. For each channel, we 



compute I K (5.2). Then, we multiply by the following LDMEs: 

(0|O4 / ' /, ( 3 5l 11 )|0) = m^v 2 ) 2 x 1.32 GeV 3 , (8.6a) 

(0|Oo /V, ( 3 P [8] )|0) = -0.109 GeV 5 , (8.6b) 

(OlOo^C^lO) = 3.12 x 10~ 3 GeV 3 , (8.6c) 

(0|Oq /V '( 1 4 81 )|0) = 4.50 x 10- 2 GeV 3 , (8.6d) 



where m 2 (v 2 ) is given in Eq. (8.5), and the numerical values on the right sides of Eq. (8.6) 
come from the fit at NLO in a s to the Tevatron and HERA data in Ref. [7]. The results 
of this computation are shown in Table 2. 

We see from the second row of Table 2 that the QQ( 3 S^) channel makes a small 
contribution to -D4 [g — > J ftp] at = 20 GeV, confirming our previous conclusion that this 
channel is not important phenomenologically at the current level of precision. We also see 
that the QQi^S^) and QQi^P^) channels give comparable contributions to D^[g — > J/ip] 
at p T = 20 GeV. As we have mentioned, at high p T at LO in a s , the fragmentation 
contribution gives the bulk of the contribution in the QQ( 3 sf^) channel (Ref. [45]). As we 
have also mentioned, the correction to the production rate in the QQ( 3 sf'^ ) channel at NLO 
in a s is only about 14% at the Tevatron [46]. Hence, we can estimate the fragmentation 
contribution to J/ip production at the Tevatron in the QQ( 3 P^) channel by multiplying 
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Table 2. Relative contributions to daJJ^ in order w 4 . The first and second rows give Io(d) for the 
short-distance coefficients in Eq. (3.f4). The second and third rows give Io{d) times the LDMEs 
in Eq. (8.6). The fourth and fifth rows give 15.2(d). The fifth and sixth rows give 15.2(d) times 
the LDMEs in Eq. (8.6). We take m = m c = 1.5 GeV. For compatibility with Ref. [7], we take 
a s = a s (m T ), where m T — \fpj> + 4mj?. We choose the point p T = 20 GeV, which implies that 
a s (m T ) = 0.154. 



I K (d)\ channel 
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10 6 GeV 3 
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10 6 LDME 
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the ratio of the QQ( 3 P^) and QQ^Sf 1 ) entries in the seventh row of Table 2 by the 
value of da/dp T x B(J/ip — > \i[i) at 20 GeV from Fig. 1(c) of Ref. [7] and by dividing by 
1.14 to account for the NLO correction in the QQ(^S^) channel. Our estimate is that the 
fragmentation contribution to da/dp? x B(J/ip — > fj,fj,) at pj- = 20 GeV from the QQ( 3 P^) 
channel is about 6 x 10 -3 nb/GeV. We see that this is comparable to (about a factor of 2 
larger than) the total NLO contribution in the QQ( 3 P^) channel in Fig. 1(c) of Ref. [7]. 
A more precise calculation of the fragmentation contribution in the QQ( 3 P^) channel will 
be necessary in order to determine whether it is the dominant contribution in that channel 
at NLO in a s at high p T . 

9 Summary 

We have calculated NRQCD short-distance coefficients for gluon fragmentation into a 3 Si 
heavy-quarkonium state through relative order v A . Our principal new result is the expres- 
sion for the sum of the relative-order-?; 4 short-distance coefficients d^\[g — > QQ( S S^)]^ + 
^4,2 [<7 ~~ ^ QQ^S^)]^ for gluon fragmentation through the 3 Si color-singlet channel. This 
expression is given in Eq. (7.7) and in the parametrization of d±\g — > QQ( 3 s[ 1 ')] fimtc in 
Eq. (B.2). As a byproduct of this calculation, we have computed the short-distance co- 
efficient for gluon fragmentation through the 3 Si color-octet channel [Eq. (7.4)], finding 
agreement with the results in Refs. [30, 31], and the sum of short-distance coefficients for 
gluon fragmentation through the 3 Pj color-octet channel [Eq. (7.5)], finding agreement 
with the result in Ref. [29] for the 3 Pj color-singlet short-distance coefficients. We have 
also computed the short-distance coefficients for gluon fragmentation through the 3 Si color- 
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singlet channel at leading order in v and at relative-order v 2 and find agreement with the 
results in Refs. [19, 20] and Ref. [21], respectively. 

The analysis in this paper involves, for the first time in an NRQCD factorization 
calculation, both single and double soft divergences in the full-QCD process. These soft 
divergences manifest themselves as single and double poles in e = (4 — d)/2 in dimen- 
sional regularization. We have dealt with the soft divergences by devising subtractions 
that remove both the single and double poles in e. We have calculated the phase-space 
integrals for the subtraction contributions analytically in d = 4 — 2e dimensions and have 
calculated the phase-space integrals for the finite remainder contributions numerically in 
d = 4 dimensions. We have also calculated the perturbative NRQCD LDMEs that appear 
in the NRQCD matching equations for the short-distance coefficients. These perturbative 
LDMEs involve both one-loop and two-loop renormalizations of the NRQCD operators. 
Our results for the renormalization-group evolution of the renormalized LDMEs confirm 
some results in Refs. [2, 23, 44], but disagree with one of the results in Ref. [23]. 

In Tables 1 and 2, we have given estimates of the relative sizes of the contributions 
of the various channels to gluon fragmentation into J/ip through relative order v 4 . As we 
have mentioned, the contribution of the 3 Si color-octet channel is believed to be important 
phenomenologically at large p T . Hence, one might expect the order-v 4 contribution in the 
3 Si color-singlet channel to be important as well, since the color-singlet channel mixes with 
the color-octet channel through single and double logarithms of the NRQCD factorization 
scale at order v 4 . Indeed, we find that the contribution to the cross section at p T = 20 GeV 
of the 3 Si color-singlet channel in order v 4 is about a factor of 2 larger than the contribution 
of the 3 Si color-singlet channel at leading order in v when the factorization scale is taken to 
be /iA = m c and is about a factor of 4 larger than the contribution of the 3 Si color-singlet 
channel at leading order in v when the factorization scale is taken to be /Ua = 2m c . This 
is a large enhancement, since one would nominally expect the order- 1; 4 contribution to be 
about (v 2 ) 2 ~ 0.04 times the order-i>° contribution. In spite of this large enhancement 
of the fragmentation contribution in order v 4 , the corresponding contribution to the J/tp 
production cross section at the Tevatron or the LHC is not important at the current level 
of precision of the phenomenology. 

We attribute the large enhancement of the order-f 4 contribution to the 3 Si color- 
singlet channel to the peaking of the sum of short-distance coefficients near z = 1. This 
peaking arises from terms that are proportional to 5(1 — z), [1/(1 — z)] + , and [log(l — 
z)/(l — z)]+. These terms are remnants of the soft divergences that appear in the full- 
QCD expression for the fragmentation process and that, ultimately, cancel in the NRQCD 
matching equations for the short-distance coefficients. We expect such a peaking, and the 
corresponding enhancement, to be present whenever soft divergences appear in a full-QCD 
process and are cancelled in the NRQCD matching equations. 7 

The signature of this peaking/enhancement in Table 2 is that the magnitude of the 

7 There is also some peaking/enhancement in the order-v 2 contribution to the 3 Si color-singlet channel. 
The order-n 2 contribution contains terms that are the product of an order-w° amplitude with an order-v 2 
amplitude. The order-u 2 amplitude is sufficiently singular that its square would produce a soft divergence. 
However, in combination with the order-«° amplitude, it produces only a peaking at z = f . 
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-Z5.2 entry is comparable to or greater than the magnitude of the Iq entry. We see this 
signature of peaking/enhancement in the 3 S*i and 3 Pj color-octet channels, as well as in 
the 3 Si color-singlet channel in order v 4 , but not in the 1 So color-octet channel. The 3 Pj 
color-octet channel contains soft divergences in the full-QCD process, but the 1 Sq color- 
octet channel does not. (The fragmentation contribution in the 3 Si color-octet channel is 
proportional to 6(1 — z) at leading order in a s .) 

The fragmentation contribution to the 3 Si color-octet channel dominates the contri- 
bution of that channel to the J/ip production cross section at hadron-hadron colliders at 
large p T , both at LO in a s (Ref. [45]) and at NLO in a s (Ref. [46]). We estimate that 
the fragmentation contribution to the 3 Pj color-octet channel gives a substantial part of 
the contribution of that channel at NLO in a s to the total J/ip production cross section 
at hadron-hadron colliders at large p T . A more precise calculation of the 3 Pj color-octet 
fragmentation contribution will be required in order to determine if it is the dominant 
contribution at large p T in that channel. 
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A Analytic calculation of phase-space integrals 

In this appendix, we give some of the details of the analytic integrations over the final-state 
phase space for £12, Si and S2. 

We first carry out the average over the angles of the transverse components of the 
final-state gluon momenta in the phase-space (4.8c). Under this angular averaging we have 
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After this angular averaging, the squared amplitude is independent of the directions of k\±_ 
and k>2±. 



In evaluating the integrals over and |&2_i_|, it is convenient to define dimensionless 
variables u\ and U2- 
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The integrals of S12, Si and S% over the magnitudes of the transverse components of 
the final-state gluon momenta can be expressed as linear combinations of the following 
elementary integrals: 
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where m takes on integer values equal to or greater than —1, and n takes on non- negative 
integer values. These integrals can be evaluated straightforwardly to obtain 
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After we have carried out the integrations over the transverse components of the final-state 
gluon momenta the remaining integrals over z\ and Z2 for S12 are simple. The integrations 
over z\ and Z2 for S\ and S2 can be expressed as linear combinations of the integrals 
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where 



/ F(z u z 2 ) = I dz! I dz 2 F( Zl ,z 2 )5{l- z- Zl - z 2 ). (A.7) 
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By making use of the identity z 2 j{z + z 2 ) = 1 — zj (z + z 2 ), we derive the relations 



D n = C n -l — zD n -\. 



(A.8a) 
(A.8b) 



Applying the recursion relations (A. 8) repeatedly, we can reduce all of the integrals to the 
forms Ar,, B„ and C„, with 1 < n < 4. It turns out that the coefficient of Dn vanishes. 



The coefficient of An is of order e°, while the coefficients of the B n are of order e _1 or e°. 
Therefore, we evaluate An only through order e°, and we evaluate the B n only through 
order e. The expressions for Aq, B n and C n can be obtained conveniently by making use 
of the identity 
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which applies when the domain of integration is < z\ < 1 — z. The distribution in the 
second term of Eq. (A.9) is defined by 
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A straightforward evaluation of the integrals then gives 
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In Eq. (A.llc), the first term in the braces is the pole contribution and remainder is finite. 
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The results for Y n for < n < 4 are 
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B Parametrizations of the short-distance coefficients 
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In this appendix, we give parametrizations for the short-distance coefficients do[g — > 
QQ( 3 5{ 1] )]( 3 )(z), d 2 [g -> QQ( 3 5[ 1] )]( 3 )(z), and d A [g -»• QQ( 3 ,s[ 1] )] finit0 (z), which are de- 
fined in Eqs. (3.11), (3.12) and (7.8), respectively. 

We observe that d [g -> QQ( 3 Si ] )]W(z), d 2 [g -»■ QQ( 3 M 1] )] (3) (^), and d 4 [p -»• QQ 
( 3 S'( 1 ^)] fimte (z) are continuous functions of 2 over the whole range < z < 1 and that they 
are not analytic at the endpoints z = and 1. One can prove that d n (0) = for n = 0, 2, 
and 4. In addition, d n (l) = a 2 b n /m 3+n is finite and calculable analytically. The values of 
b n for n = 0, 2, and 4 are 
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Using this information, we parametrize do[g ->■ QQ( 3 5'{ 1 ')]( 3 )(z), a! 2 [g ->■ QQ^S 1 ^ 1 )]^ (z), 
and d 4 [ 5 -> QQ( 3 sf ] )f nitc (z) as follows: 
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In order to fix the parameters in Eq. (B.2), we have computed dob ~~ y QQ( 3 ^i^)]^ 3 K z )i 
d 2 [g -> QQ{ 3 Si ] )]W(z), and d 4 [g -> QQ( 3 sI 1] )] finitc (z) by numerical integration for 100 
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Table 3. Results of fitting the parametrization in Eq. (B.2) to numerical values for the color- 
singlet short-distance coefficients d [g -t QQ( 3 Si^)]^(z), d 2 [g -> QQ( 3 5'[ 1 ')] (3) (z), and d A [g -> 
QQ{ 3 S l i ] )f nitc {z), which are defined in Eqs. (3.11), (3.12) and (7.8), respectively. The quantity 
AJL* is described in the text. 
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points in each of the ranges z = to z = 10 -2 , z = 10 -2 to z = 1 — 10~ 2 , and z = 1 — 10~ 2 
to z = 1. We have carried out the integrations in two ways: (1) using the change of 
variables that was proposed in Refs. [19, 20] and is described in Appendix B of Ref. [21] 
and (2) using the change of variables that is given in Eq. (A. 2). At each value of z, the 
numerical results from the two methods of integration agree to better than one part in 10 3 . 

We fit the parametrizations to these points, using x 2 = J2i[d( z i) ~d fit (zi)] 2 /[a(zi)] 2 , as 
the criterion for goodness of fit, where a{zi) is the error in the numerical integration that 
is given by the VEGAS integration program [48]. In cases in which it is possible to reduce 
the number of parameters in the fit, without significantly affecting x 2 , we have done so. 

In Table 3, we show the results of this fitting procedure, along with the value for each 
fit of A™^ = M&x\d(zi) — d fit (zj)|. In these fits, \ 2 P er degree of freedom is much larger 
than one because the error in the fit is generally considerably larger than the errors in the 
numerical integrations that are given by VEGAS. A conservative estimate of the overall 
error in the parametrization at each value of z is y^A^J 2 + [10~ 3 x d n (z)] 2 . 
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